1 C Q-D PROGRAM
2 DIMENSION NDLTH(4),CDH(4),ARR(5,5),AOR(5,5),ARO(5,5),AOO(5,5),
3 1BRR(5,5),BOR(5,5),BRO(5,5),BOO(5,5),CRR(5,5),COR(5,5),CRO(5,5),
4 1COO(5,5),DRR(5,5),DOR(5,5),DRO(5,5),DOO(5,5),PN(7),X1(6),X2(6),
5 1Y1(6),Y2(6),Q(10,10)
6 DOUBLE PRECISION THETA,CDH,DLTH,SRMUL,STH,CTH,PN,DPNI,DPNJ,FCT,B1,
7 1B2,C1,C2,PMY1,PMY2,GA,R,DR,NIND,KAPPA,ALPHA,A,B,C,D,PI,THETA1,
8 1THETA2,THETA3,THETA4,E,F,G
9 COMPLEX*16 CI,K1R,K2R,FMY,DK1R,DK2R,X1,X2,Y1,Y2,Z1I,Z2J,DKRX1I,
10 1DKRY1I,DKRZ1I,DKRX2J,DKRY2J,DKRZ2J,W1,W2,W3,W4,W5,W6,W7,W8,W9,
11 1ARR,AOR,ARO,AOO,BRR,BOR,BRO,BOO,CRR,COR,CRO,COO,DRR,DOR,DRO,DOO,
12 1K1,K2,Q,FC
13 COMMON/FAC/FCT (100),NRISK,NTRIX
14 NAMELIST/HEJ/ARR,AOR,ARO,AOO,BRR,BOR,BRO,BOO,CRR,COR,CRO,COO,DRR,
15 1DOR,DRO,DOO/NEJ/Q
16 cc CALL UNSTAE
17 PI = 4.0D0*DATAN(1.0D0)
18 NRISK = 57
19 NTRIX = 39
20 NEND = 78
21
22 c half the diamention of Q -matrix NR > 2|k|a
23 NR = 5
24
25 c polarization index 0=orthogonal, 1=parallel
26 NP = 0
27
28 c paramter 0=homogeneous, 1=layered
29 LAY = 0
30
31 c permeability outside
32 PMY1 = 1.0D0
33
34 c permeability inside
35 PMY2 = 1.0D0
36
37 c complex wave vektor outside
38 K1 = (0.75D0,0.0D0)
39
40 c parameters for wave vektor inside
41 c in the way k2= nind*(1+i*kappa)*k1, Im(k2) > 0
42 NIND = 4.0D0/5.0D0
43 KAPPA = 0.0D0
44 C = -0.1D0
45 c size parameter
46 A = 1.0D0
47
48 c number of sections in the integration
49 NSECT = 1
50
51 c number of integration intervalls in section I
52 c dividible by 4
53 NDLTH(1) = 192
54
55 c endpoint in section I resp. starting point in section I+1
56 c of integration (rad)
57 THETA1 = PI
58
59 CDH(1) = THETA1/DFLOAT(NDLTH(1))
60 N = NRISK
61 L = NTRIX
62 M = NEND-2
63 N1 = N-1
64 DO 5 K = 1,100
65 5 FCT(K) = 0.0D0
66 FCT(1) = 1.0D0
67 DO 6 K = 1,N1
68 6 FCT(K+1) = FCT(K)*DFLOAT(K)
69 FCT(N+1) = 0.1D0**L*FCT(N)*DFLOAT(N)
70 DO 7 K = N,M
71 7 FCT (K+2) = FCT(K+1)*DFLOAT(K+1)
72 FC = (1.0D0,0.0D0)
73 IF(NP.EQ.1) FC = -FC
74 CI = (0.0D0, 1.000)
75 FMY = DCMPLX(PMY1/PMY2,0.0D0)
76 K2 = DCMPLX(NIND,NIND*KAPPA)*K1
77 ND = 2*NR
78 NR1 = NR+1
79 NR2 = NR+2
80
81 iq=0
82 DO 150 M1 = 1,NR1
83 M = M1-1
84 MD = M
85 IF(M.EQ.0) MD = 1
86 DO 25 J = 1,NR
87 DO 25 I = 1,NR
88 C Seite 5
89
90
91 ARR(I,J) = (0.0D0,0.0D0)
92 AOR(I,J) = (0.0D0,0.0D0)
93 ARO(I,J) = (0.0D0,0.0D0)
94 AOO(I,J) = (0.0D0,0.0D0)
95 BRR(I,J) = (0.0D0,0.0D0)
96 BOR(I,J) = (0.0D0,0.0D0)
97 BRO(I,J) = (0.0D0,0.0D0)
98 BOO(I,J) = (0.0D0,0.0D0)
99 CRR(I,J) = (0.0D0,0.0D0)
100 COR(I,J) = (0.0D0,0.0D0)
101 CRO(I,J) = (0.0D0,0.0D0)
102 COO(I,J) = (0.0D0,0.0D0)
103 DRR(I,J) = (0.0D0,0.0D0)
104 DOR(I,J) = (0.0D0,0.0D0)
105 DRO(I,J) = (0.0D0,0.0D0)
106 DOO(I,J) = (0.0D0,0.0D0)
107 25 CONTINUE
108 THETA = 0.0D0
109 DO 90 ISECT = 1,NSECT
110 NTHETA = NDLTH(ISECT)+1
111 DLTH = CDH(ISECT)
112 NFIX = 1
113 DO 89 K = 1,NTHETA
114 IF(K-1)31,31,32
115 31 CONTINUE
116 SRMUL = DLTH*7.0D0/22.5D0
117 IF(ISECT.EQ.1) GO TO 89
118 GO TO 41
119 32 CONTINUE
120 IF(K-NTHETA)34,33,33
121 33 CONTINUE
122 SRMUL = DLTH*7.0D0/22.5D0
123 IF(ISECT-NSECT)40,89,89
124 34 CONTINUE
125 GO TO (35,36,37,38),NFIX
126 35 CONTINUE
127 SRMUL = DLTH*32.0D0/22.5D0
128 NFIX = 2
129 GO TO 39
130 36 CONTINUE
131 SRMUL = DLTH*12.0D0/22.5D0
132 NFIX = 3
133 GO TO 39
134 37 CONTINUE
135 SRMUL = DLTH*32.0D0/22.5D0
136 NFIX = 4
137 GO TO 39
138 38 CONTINUE
139 SRMUL = DLTH*14.0D0/22.5D0
140 NFIX = 1
141 39 CONTINUE
142 40 CONTINUE
143 THETA = THETA+DLTH
144 STH = DSIN(THETA)
145 CTH = DCOS(THETA)
146 CALL LEG(THETA,M,NR2,PN)
147 C Seite 6
148
149
150 41 CONTINUE
151 CALL TRCIR(STH,CTH,C,A,R,DR)
152 K1R= K1*DCMPLX(R,0.0D0)
153 DK1R = K1*DCMPLX(DR,0.0D0)
154
155
156 K2R = K2*DCMPLX(R,0.0D0)
157 IF(K.EQ.1) GO TO 52
158 CALL CBN(K1R,NR1,X1,Y1)
159 CALL CBN(K2R,NR1,X2,Y2)
160 B1 = CDABS(K1R*K1R*(X1(2)*Y1(1)-X1(1)*Y1(2))-DCMPLX(1.0D0,0.0D0))
161 B2 = CDABS(K2R*K2R*(X2(2)*Y2(1)-X2(1)*Y2(2))-DCMPLX(1.0D0,0.0D0))
162 C1 = CDABS(K1R*K1R*(X1(NR1)*Y1(NR)-X1(NR)*Y1(NR1))-DCMPLX(1.0D0,0.
163 10D0))
164 C2 = CDABS(K2R*K2R*(X2(NR1)*Y2(NR)-X2(NR)*Y2(NR1))-DCMPLX(1.0D0,0.
165 10D0))
166 IF(B1.LT.1.D-10.AND.C1.LT.1.D-10) GO TO 47
167 PRINT 45
168 45 FORMAT('0 ERROR IN BESSEL-NEUMAN-1 TEST')
169 PRINT 46, B1,C1,ISECT,K
170 46 FORMAT(2D25.15,2I5)
171 47 CONTINUE
172 IF(B2.LT.1.D-10.AND.C2.LT.1.D-10) GO TO 52
173 PRINT 48
174 48 FORMAT ('0 ERROR IN BESSEL-NEUMAN-2 TEST')
175 49 FORMAT(2D25.15,2I5)
176 PRINT 49,B2,C2,ISECT,K
177 52 CONTINUE
178 DO 88 I = MD,NR
179 DO 88 J = MD,NR
180 I1 = I+J
181 J1 = J+1
182 I2 = I + 2
183 J2 = J+2
184 DPNI = -(DFLOAT(I1)*CTH*PN(I1)-DFLOAT(I1-M)*PN(I2))/STH
185 DPNJ = -(DFLOAT(J1)*CTH*PN(J1)-DFLOAT(J1-M)*PN(J2))/STH
186 Z1I = X1(I1)+CI*Y1(I1)
187 DKRX1I = K1R*X1(I)-DCMPLX(DFLOAT(I),0.0D0)*X1(I1)
188 DKRY1I = K1R*Y1(I)-DCMPLX(DFLOAT(I),0.0D0)*Y1(I1)
189 DKRZ1I = DKRX1I+CI*DKRY1I
190 Z2J = X2(J1)*CI*Y2(J1)
191 DKRX2J = K2R*X2(J)-DCMPLX(DFLOAT(J),0.0D0)*X2(J1)
192 DKRY2J = K2R*Y2(J)-DCMPLX(DFLOAT(J),0.0D0)*Y2(J1)
193 DKRZ2J = DKRX2J+CI*DKRY2J
194 W1 = DCMPLX(SRMUL*STH*(DPNI*DPNJ*DFLOAT(M*M)*PN(I1)*PN(J1)/STH**2)
195 1,0.0D0)
196 W2 = DCMPLX(SRMUL*STH*DFLOAT(I*I1)*PN(I1)*DPNJ,0.0D0)
197 W3 = DCMPLX(SRMUL*STH*DFLOAT(J*J1)*PN(J1)*DPNI,0.0D0)
198 ARR(I,J) = ARR(I,J)+K1R*DKRX1I*X2(J1)*W1+DK1R*X1(I1)*X2(J1)*W2
199 AOR(I,J) = AOR(I,J)+K1R*DKRZ1I*X2(J1)*W1+DK1R*Z1I*X2(J1)*W2
200 DRR(I,J) = DRR(I,J)+K1R*X1(I1)*DKRX2J*W1+DK1R*X1(I1)*X2(J1)*W3
201 DOR(I,J) = DOR(I,J)+K1R*Z1I*DKRX2J*W1+DK1R*Z1I*X2(J1)*W3
202 IF(LAY.EQ.0) GO TO 53
203 ARO(I,J) = ARO(I,J)+K1R*DKRX1I*Z2J*W1+DK1R*X1(I1)*Z2J*W2
204 AOO(I,J) = AOO(I,J)+K1R*DKRZ1I*Z2J*W1*DK1R*Z1I*Z2J*W2
205 DRO(I,J) = DRO(I,J)+K1R*X1(I1)*DKRZ2J*W1+DK1R*X1(I1)*Z2J*W3
206 DOO(I,J) = DOO(I,J)+K1R*Z1I*DKRZ2J*W1+DK1R*Z1I*Z2J*W3
207 53 CONTINUE
208 IF(M.EQ.0) GO TO 88
209 C Seite 7
210
211
212 W4 = DCMPLX(SRMUL*(DPNI*PN(J1)+PN(I1)*DPNJ),0.0D0)
213 W5 = DCMPLX(SRMUL*PN(I1)*PN(J1),0.0D0)
214 BRR(I,J) = BRR(I,J)+K1R*K1R*X1(I1)*X2(J1)*W4
215 BOR(I,J) = BOR(I,J)+K1R*K1R*Z1I*X2(J1)*W4
216 W6 = DK1R*(X1(I1)*DKRX2J*DCMPLX(DFLOAT(I*(I+1)),0.0D0)+DCMPLX(
217 1DFLOAT(J*(J+1)),0.0D0)*DKRX1I*X2(J1))/K1R
218 CRR(I,J) = CRR(I,J)+DKRX1I*DKRX2J*W4+W5*W6
219 W6 = DK1R*(Z1I*DKRX2J*DCMPLX(DFLOAT (I*(I+1)),0.0D0)+DCMPLX(
220 1DFLOAT(J*(J+1)),0.0D0)*DKRZ1I*X2(J1))/K1R
221 COR(I,J) = COR(I,J)+DKRZ1I*DKRX2J*W4+W5*W6
222 IF(LAY.EQ.0) GO TO 88
223 BRO(I,J) = BRO(I,J)+K1R*K1R*X1(I1)*Z2J*W4
224 BOO(I,J) = BOO(I,J)+K1R*K1R*Z1I*Z2J*W4
225 W6 = DK1R*(X1(I1)*DKRZ2J*DCMPLX(DFLOAT(I*(I+1)),0.0D0)+DCMPLX(
226 1DFLOAT(J*(J+1)),0.0D0)*DKRX1I*Z2J)/K1R
227 CRO(I,J) = CRO(I,J)+DKRX1I*DKRZ2J*W4+W5*W6
228 W6 = DK1R*(Z1I*DKRZ2J*DCMPLX(DFLOAT(I*(I+1)),0.0D0)+DCMPLX(
229 1DFLOAT(J*(J+1)),0.0D0)*DKRZ1I*Z2J)/K1R
230 COO(I,J) = COO(I,J)+DKRZ1I*DKRZ2J*W4+W5*W6
231 88 CONTINUE
232 89 CONTINUE
233 90 CONTINUE
234 DO 98 I = MD,NR
235 DO 98 J = MD,NR
236 I1 = I+1
237 J1 = J+1
238 GA = DSQRT(DFLOAT((2*I+1)*(2*J+1))/DFLOAT(I*I1*J*J1))/2.0D0
239 IF(M)96,96,95
240 95 CONTINUE
241 GA = GA*DSQRT(FCT(I1-M)*FCT(J1-M)/FCT(I1+M)*FCT(J1*M))
242 96 CONTINUE
243 W1 = DCMPLX(GA,0.0D0)
244 ARR(I,J) = W1*ARR(I,J)
245 AOR(I,J) = W1*AOR(I,J)
246 DRR(I,J) = W1*DRR(I,J)
247 DOR(I,J) = W1*DOR(I,J)
248 IF(LAY.EQ.0) GO TO 97
249 ARO(I,J) = W1*ARO(I,J)
250 AOO(I,J) = W1*AOO(I,J)
251 DRO(I,J) = W1*DRO(I,J)
252 DOO(I,J) = W1*DOO(I,J)
253 97 CONTINUE
254 IF(M.EQ.0) GO TO 98
255 W1 = DCMPLX(DFLOAT(M),0.0D0)*W1
256 BRR(I,J) = W1*BRR(I,J)
257 BOR(I,J) = W1*BOR(I,J)
258 CRR(I,J) = W1*CRR(I,J)
259 COR(I,J) = W1*COR(I,J)
260 IF(LAY.EQ.0) GO TO 98
261 BRO(I,J) = W1*BRO(I,J)
262 BOO(I,J) = W1*BOO(I,J)
263 CRO(I,J) = W1*CRO(I,J)
264 COO(I,J) = W1*COO(I,J)
265 98 CONTINUE
266 IF(M.GT.3) GO TO 99
267 WRITE(6,HEJ)
268 99 CONTINUE
269 C Seite 8
270
271
272 LLAY = 2+LAY*2
273 DO 120 LL = 1,LLAY
274 DO 110 I = 1,NR
275 DO 110 J = 1,NR
276 GO TO (101,102,103,104), LL
277 101 CONTINUE
278 Q(2*I-1,2*J-1) = -ARR(I,J)+FMY*DRR(I,J)
279 Q(2*I-1,2*J) = FC*(K1*CRR(I,J)/K2+K2*FMY*BRR(I,J)/K1)
280 Q(2*I,2*J-1) = -FC*(BRR(I,J)+FMY*CRR(I,J))
281 Q(2*I,2*J) = (K1*DRR(I,J)/K2-K2*FMY*ARR(I,J)/K1)
282
283
284 GO TO 110
285 102 CONTINUE
286 Q(2*I-1,2*J-1) =-AOR(I,J)+FMY*DOR(I,J)
287 Q(2*I-1,2*J) = FC*K1*COR(I,J)/K2+K2*FMY*BOR(I,J)/K1
288 Q(2*I,2*J-1) = -FC*(BOR(I,J)*FMY*COR(I,J))
289 Q(2*I,2*J) = (K1*DOR(I,J)/K2-K2*FMY*AOR(I,J)/K1)
290 GO TO 110
291 103 CONTINUE
292 Q(2*I-1,2*J-1) = -ARO(I,J)+FMY*DRO(I,J)
293 Q(2*I-1,2*J) = FC*(K1*CRO(I,J)/K2+K2*FMY*BRO(I,J)/K1)
294 Q(2*I,2*J-1) = -FC*(BRO(I,J)+FMY*CRO(I,J))
295 Q(2*I,2*J) = (K1*DRO(I,J)/K2-K2*FMY*AOO(I,J)/K1)
296 104 CONTINUE
297 Q(2*I-1,2*J-1) = -AOO(I,J)+FMY*DOO(I,J)
298 Q(2*I-1,2*J) = FC*(K1*COO(I,J)/K2+K2*FMY*BOO(I,J)/K1)
299 Q(2*I,2*J-1) = -FC*(BOO(I,J)+FMY*COO(I,J))
300 Q(2*I,2*J) = (K1*DOO(I,J)/K2-K2*FMY*AOO(I,J)/K1)
301 110 CONTINUE
302 WRITE (6,NEJ)
303 GO TO (111,111,112,112)LL
304 111 CONTINUE
305
306 iq=iq+1
307 WRITE (21+iq) Q
308 ENDFILE 21+iq
309 GO TO 120
310 112 CONTINUE
311 WRITE(32) Q
312 ENDFILE 32
313 120 CONTINUE
314 PRINT 777, M1
315 777 FORMAT(I4)
316 150 CONTINUE
317 PRINT 199
318 199 FORMAT('0 Q-MATRICES NOW HOPEFULLY PLACED INTO TAPE')
319 STOP
320 cc DEBUG SUBCHK
321 END
322 C Seite 9
323
324 SUBROUTINE BESSEL(NORDER,ARGMNT,ANSWR,IERROR)
325 DOUBLE PRECISION ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
326 IERROR = 0
327 N = NORDER
328 X = ARGMNT
329 SUM = 1.0D0
330 APR = 1.0D0
331 TOPR = -0.5D0*X*X
332 CI = 1.0D0
333 CNI = DFLOAT(2*N+3)
334 DO 60 I = 1,100
335 ACR = TOPR*APR/(CI*CNI)
336 SUM = SUM+ACR
337 IF(DABS(ACR/SUM)-1.0D-20)100,100,40
338 40 APR = ACR
339 CI = CI+1.0D0
340 CNI = CNI+2.0D0
341 60 CONTINUE
342 IERROR = I
343 PRINT 10
344 10 FORMAT(24HO ERROR IN SUM OF BESSEL)
345 GO TO 200
346 100 PROD = DFLOAT(2*N+1)
347 FACT = 1.0D0
348 IF(N)160,160,120
349 120 DO 140 IFCT = 1,N
350 FACT = FACT*X/PROD
351 PROD = PROD-2.0D0
352 140 CONTINUE
353 160 ANSWR = FACT*SUM
354 200 RETURN
355 END
356 C Seite 23
357
358
359 SUBROUTINE BN(PCKR,NRINK,BSSLSP,CNEUMN)
360 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
361 DOUBLE PRECISION PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,
362 1SNSC,BSSLSP,CNEUMN
363 NRANKI = NRINK
364 NRANK = NRANKI-1
365 NVAL = NRANK-1
366 DO 40 I = 1,4
367 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
368 IF(IERROR)20,20,32
369 20 ANSA = ANSWR
370 NVAL = NVAL+1
371 CALL BESSEL(NVAL,PCKR,ANSWR,IERROR)
372 IF(IERROR)24,24,28
373 24 ANSB = ANSWR
374 GO TO 60
375 28 NVAL = NVAL-1
376 32 NVAL = NVAL+NRANK
377 40 CONTINUE
378 STOP
379 60 IF(NVAL-NRANK)100,100,64
380 64 IEND = NVAL-NRANK
381 CONN = DFLOAT(2*(NVAL-1)+1)
382 DO 72 IP = 1,IEND
383 ANSC = CONN*ANSA/PCKR-ANSB
384 CONN = CONN-2.0D0
385 ANSB = ANSA
386 ANSA = ANSC
387 72 CONTINUE
388 100 BSSLSP(NRANKI) = ANSB
389 BSSLSP(NRANKI-1) = ANSA
390 CONN = DFLOAT(NRANK+NRANK-1)
391 IE = NRANKI-2
392 JE = IE
393 DO 120 JB = 1,JE
394 ANSC = CONN*ANSA/PCKR-ANSB
395 BSSLSP(IE) = ANSC
396 ANSB = ANSA
397 ANSA = ANSC
398 IE = IE-1
399 CONN = CONN-2.0D0
400 120 CONTINUE
401 CMULN = 3.0D0
402 SNSA = -DCOS(PCKR)/PCKR
403 SNSB = -DCOS(PCKR)/(PCKR*PCKR)-DSIN(PCKR)/PCKR
404 CNEUMN(1) = SNSA
405 CNEUMN(2) = SNSB
406 DO 280 I = 3,NRANKI
407 SNSC = CMULN*SNSB/PCKR-SNSA
408 CNEUMN(I) = SNSC
409 SNSA = SNSB
410 SNSB = SNSC
411 CMULN = CMULN+2.0D0
412 280 CONTINUE
413 RETURN
414 END
415 C Seite 24
416
417
418 SUBROUTINE CBESS (NORDER,ARGMNT,ANSWR,IERROR)
419 COMPLEX*16 ARGMNT,ANSWR,X,SUM,APR,TOPR,CI,CNI,ACR,PROD,FACT
420 IERROR = 0
421 N = NORDER
422 X = ARGMNT
423 SUM = (1.0D0,0.0D0)
424 APR = (1.0D0,0.0D0)
425 TOPR = -(0.5D0,0.0D0)*X*X
426 CI = (1.0D0,0.0D0)
427 CNI = DCMPLX(DFLOAT(2*N+3),0.0D0)
428 DO 60 I = 1,100
429 ACR = TOPR*APR/(CI*CNI)
430 SUM = SUM+ACR
431 IF(CDABS(ACR/SUM)-1.0D-20)100,100,40
432 40 APR = ACR
433 CI = CI+(1.0D0,0.0D0)
434 CNI = CNI+(2.0D0,0.0D0)
435 60 CONTINUE
436 IERROR = 1
437 PRINT 10
438 10 FORMAT('0 ERROR IN SUM OF CBESS')
439 GO TO 200
440 100 PROD = DCMPLX(DFLOAT(2*N+1),0.0D0)
441 FACT = (1.0D0,0.0D0)
442 IF(N)160,160,120
443 120 DO 140 IFCT = 1,N
444 FACT = FACT*X/PROD
445 PROD = PROD-(2.0D0,0.0D0)
446 140 CONTINUE
447 160 ANSWR = FACT*SUM
448 200 RETURN
449 END
450 C Seite 25
451
452
453 SUBROUTINE CBN(PCKR,NRINK,BSSLSP,CNEUMN)
454 DIMENSION BSSLSP(NRINK),CNEUMN(NRINK)
455 COMPLEX*16 PCKR,ANSWR,ANSA,ANSB,ANSC,CONN,CMULN,SNSA,SNSB,SNSC,
456 1BSSLSP,CNEUMN
457 NRANKI = NRINK
458 NRANK = NRANKI-1
459 NVAL = NRANK-1
460 DO 40 I = 1,4
461 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
462 IF(IERROR)20,20,32
463 20 ANSA = ANSWR
464 NVAL = NVAL+1
465 CALL CBESS(NVAL,PCKR,ANSWR,IERROR)
466 IF(IERROR)24,24,28
467 24 ANSB = ANSWR
468 GO TO 60
469 28 NVAL = NVAL-1
470 32 NVAL = NVAL+NRANK
471 40 CONTINUE
472 STOP
473 60 IF(NVAL-NRANK)100,100,64
474 64 IEND = NVAL-NRANK
475 CONN = DCMPLX(DFLOAT(2*(NVAL-1)+1),0.0D0)
476 DO 72 IP = 1,IEND
477 ANSC = CONN*ANSA/PCKR-ANSB
478 CONN = CONN-(2.0D0,0.0D0)
479 ANSB = ANSA
480 ANSA = ANSC
481 72 CONTINUE
482 100 BSSLSP(NRANKI) = ANSB
483 BSSLSP(NRANKI-1) = ANSA
484 CONN = DCMPLX(DFLOAT(NRANK+NRANK-1),0.0D0)
485 IE = NRANKI-2
486 JE = IE
487 DO 120 JB = 1,JE
488 ANSC = CONN*ANSA/PCKR-ANSB
489 BSSLSP(IE) = ANSC
490 ANSB = ANSA
491 ANSA = ANSC
492 IE = IE-1
493 CONN = CONN-(2.0D0,0.0D0)
494 120 CONTINUE
495 CMULN = (3.0D0,0.0D0)
496 SNSA = -CDCOS(PCKR)/PCKR
497 SNSB = -CDCOS(PCKR)/(PCKR*PCKR)-CDSIN(PCKR)/PCKR
498 CNEUMN(1) = SNSA
499 CNEUMN(2) = SNSB
500 DO 280 I = 3,NRANKI
501 SNSC = CMULN*SNSB/PCKR-SNSA
502 CNEUMN(I) = SNSC
503 SNSA = SNSB
504 SNSB = SNSC
505 CMULN = CMULN+(2.0D0,0.0D0)
506 280 CONTINUE
507 RETURN
508 END
509 C Seite 26
510
511
512 SUBROUTINE LEG(THETA,M,NRJNK,PNMLLG)
513 DIMENSION PNMLLG(NRJNK)
514 DOUBLE PRECISION THETA,PLA,PLB,PLC,DTWM,CNM,CNMM,CNMUL,PRODM,
515 1QUANM,PNMLLG
516 NRANKI = NRJNK
517 KMV = M
518 DTWM = DFLOAT (2*M+1)
519 IF(THETA)16,4,16
520 4 IF(KMV)12,12,6
521 6 DO 8 ILG = 1,NRANKI
522 PNMLLG(ILG) = 0.0D0
523 8 CONTINUE
524 GO TO 88
525 12 PNMLLG(1) = 1.0D0
526 PLA = 1.0D0
527 GO TO 48
528 16 IF(KMV)20,20,40
529 20 PLA = 1.0D0
530 PLB = DCOS(THETA)*PLA
531 PNMLLG(1) = PLA
532 PNMLLG(2) = PLB
533 IBEG = 3
534 GO TO 60
535 40 DO 44 ILG = 1,KMV
536 PNMLLG(ILG) = 0.0D0
537 44 CONTINUE
538 PRODM = 1.0D0
539 QUANM = DFLOAT(KMV)
540 DO 52 IFCT = 1,KMV
541 QUANM = QUANM+1.0D0
542 PRODM = QUANM*PRODM/2.0D0
543 52 CONTINUE
544 PLA = PRODM*DSIN(THETA)**KMV
545 PNMLLG(KMV+1) = PLA
546 48 PLB = DTWM*DCOS(THETA)*PLA
547 PNMLLG(KMV+2) = PLB
548 IBEG = KMV+3
549 60 CNMUL = DFLOAT(IBEG+IBEG-3)
550 CNM = 2.0D0
551 CNMM = DTWM
552 DO 80 ILGR = IBEG,NRANKI
553 PLC = (CNMUL*DCOS(THETA)*PLB-CNMM*PLA)/CNM
554 PNMLLG(ILGR) = PLC
555 PLA = PLB
556 PLB = PLC
557 CNMUL = CNMUL+2.0D0
558 CNM = CNM+1.0D0
559 CNMM = CNMM+1.0D0
560 80 CONTINUE
561 88 RETURN
562 END
563 C Seite 27
564
565 SUBROUTINE TRCIR(STH,CTH,C,A,R,DR)
566 DOUBLE PRECISION STH,CTH,C,A,R,DR,X,ST,CT
567 ST = C*STH
568 CT = C*CTH
569 X = DSQRT(A**2-ST**2)
570 R = CT+X
571 DR = -ST-ST*CT/X
572 RETURN
573 END
574
575
576 SUBROUTINE ELLIPS(STH,CTH,AZ,BRA,R,DR)
577 DOUBLE PRECISION STH,CTH,AZ,BRA,R,DR
578 R = AZ/DSQRT (CTH**2+(AZ*STH/BRA)**2)
579 DR = R**3*STH*CTH*(1.0D0-(AZ/BRA)**2)/AZ**2
580 RETURN
581 END
582
583
584 SUBROUTINE TRELLI(STH,CTH,AZ,BRA,C,R,DR)
585 DOUBLE PRECISION STH,CTH,AZ,BRA,C,R,DR,NE,RO
586 NE = (AZ*STH)**2+(BRA*CTH)**2
587 RO = NE-(C*STH)**2
588 RO = DSQRT(RO)
589 R = (BRA**2*C*CTH+AZ*BRA*RO)/NE
590 DR = -2.0D0*STH*CTH*(AZ**2-BRA**2)*R/NE+(-BRA**2*C*STH+STH*CTH*AZ*
591 1BRA*(AZ**2-BRA**2-C**2)/RO)/NE
592 RETURN
593 END
594 C Seite 40
595